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ABSTRACT 





o 


The cooling phase of thermonuclear (type-I) X-ray bursts can be used to constrain the neutron star (NS) compactness by comparing 
the observed cooling tracks of bursts to accurate theoretical atmosphere model calculations. By applying the so-called cooling tail 
method, where the information from the whole cooling track is used, we constrain the mass, radius, and distance for three different 
NSs in low-mass X-ray binaries 4U 1702-429, 4U 1724-307, and SAX J1810.8-260. Care is taken to only use the hard state bursts 
where it is thought that only the NS surface alone is emitting. We then utilize a Markov chain Monte Carlo algorithm within a Bayesian 
framework to obtain a parameterized equation of state (EoS) of cold dense matter from our initial mass and radius constraints. This 
allows us to set limits on various nuclear parameters and to constrain an empirical pressure-density relation for the dense matter. Our 
predicted EoS results in NS radius between 10.5 - 12.8 km (95% confidence limits) for a mass of 1.4 M© depending slightly on the 
assumed composition. Due to systematic errors and uncertainty in the composition these results should be interpreted as lower limits 
for the radius. 
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The equation of state (EoS) of the cold dense matter inside 
neutron stars (NS) has remained a mystery for decades. Exper¬ 
iments on Earth and theoretical many-body calculations have 
constrained the pressure-density relation of matter near the nu¬ 
clear saturation densities. Recently, progress has also been made 


in measuring the NS radii (for a review, see Miller||2013[ 


Ozel 


X 


|2013| [Suleimanov et al.|2015| ) that allow us to constrain the be¬ 
havior of the EoS at higher densities by inverting the Tolman- 
Volkoff-Qppe nheimer (TOV; |Tolman| ri939[ [Oppenheimer &| 
|Volkoff|1939| ) structure equations. Eurthermore, these measure¬ 
ments probe the phase diagram of dense quantum chromody¬ 
namics at lower temperatures and higher baryon densities than 
the measurements of, for example, ultrarela tivistic heavy-i on 
collisions inside earthly laboratories (see, e.g., |Lattimer|2Q12| ). 


One of the most promising candidates for obtaining accu¬ 
rate astrophysical mass-radius (M - R) measurements has been 
the thermal emission originating from NS surface layers. One 
possibility is to use the cooling of NS surface during type-I X- 
ray bursts from low-mass X-ray binary (LMXB) systems, where 
the cooling tail is shown to follow theoretical model predictions 
surprisingly well ( [Poutanen et al.|[20T4) . In these systems, the 
NS is accompanied by a lighter, usually a main sequence or 
evolved late-type, star that fills its Roche lobe and transfers ma¬ 
terial through accretion disk onto the NS. After accumulating 
enough material, the fuel is rapidly burned in a thermonuclear 


explosion occurring below the surface in the NS ocean. Some 
of these bursts can be so energetic that the Eddington limit is 
reached, causing the entire NS photosphere to expand. These 
photospheric radius expansion (PRE) bursts can then be used to 
obtain mass and radius measurements by comparing the cool¬ 
ing tail of the b urst to accurate theoretical predictions (for early 
work, see, e.g., Damen et al.||1990] [van Paradijs et al.]|1990t 
|Lewin et al.|1993|). 

R ecent studies ([Suleimanov et al.|[^lla[ [Poutanen et al.| 
|2014t[Kajava et al.|2014| ) have demonstrated that the X-ray burst 
cooling properties heavily depend on the accretion rate and spec¬ 
tral state of the source. The key finding was that care must be 
taken to select only those bursts that show “passive cooling", i.e. 
the ones that occur at the hard spectral state and at small accre¬ 
tion rate, where the extra heating from the in-falling material ap¬ 
pears to be negligible. [Kajava et al.| ( |2014] ) also showed that the 
evolution of the blackbody normalization can be used as a trace 
to pin down the passively cooling bursts used in the M - R mea¬ 
surements. The soft state bursts, on the other hand, show only 
weak or completely non-existent evolution of the normalization 
that is in contradiction with the theoretical atmosphere model 
predictions. 


In addition to only using the passively cooling bursts, it is 
possible to improve the analysis by using the information from 
the whole cooling tail by applying the so-called “cooling tail 
method”. In this recently developed method, the observed cool¬ 
ing track in the blackbody normalization K vs. the flux F (or 
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rather vs. F) plane is compared to the theoretical model 
evolution of the color-correction factor versus the luminosity (in 
units of the Eddington), /c - L/LEdd. that is the so-called color- 
correction curve. By comparing the whole cooling track to the 
models (with variable /c) - in contrast for example to the “touch¬ 
down method” where /c is assumed to be constant - we can infer 
more robust constraints from the data. This also allows us to cir¬ 
cumvent the problematic issue of deciding where the Eddington 
limit is reached and when the photospheric radius coincides with 
the NS radius (see, e.g., [Steiner et al.||2010| ). In fact, because it 
is the curvature of the evolving color-correction factor that is 
used to constrain the Eddington flux, the method is val id even 
for bursts that do not reach the Eddington limit (see, e.g., |Zamfir| 
et al.|2Q12| ). So far, no work exists where the cooling tail method 
would have been applied with this kind of strict hard-state burst 
selection criteria for many sources simultaneously. Therefore, 
we now pay special attention to the burst selection and choose 
only the most well-behaved hard-state PRE bursts for our anal¬ 
ysis. We then use these bursts to constrain the mass and radii of 
three different NSs by applying the cooling tail method to them. 

Using these constraints we can then also go one step further 
and address the issue of unknown EoS of the cold dense mat¬ 
ter inside neutron stars. To do this, we use Bayesian inference 
to derive empirical pressure-density and mass-radius relations 
based on our burst results. Using a Markov Chain Monte Carlo 
(MCMC) algorithm, we jointly fit the three M - R constraints 
from each source allowing us to put astrophysical constraints on 
some of the nuclear physics parameters, such as the symmetry 
energy S and t he pressure of neutron-ric h matter at the satura¬ 
tion density X ( jLattimer & Steiner|2014 ). In addition, the com¬ 
bination of the cooling tail observations and parameterized EoS 
allows us to make more accurate mass and radius measurements 
for each of the sources, indicating a new way of probing individ¬ 
ual NS characteristics. 

The paper is structured as follows. In Section]^ we present 
the methods used for the data reduction of the bursts. Then, in 
Section we use this data to obtain separate mass, radius, and 
distance constraints for the three sources in our sample. In the 
second part of the paper, in Section]^ we use Bayesian analysis 
to obtain the parameterized EoS. Einally, in the last Sectionj^ we 
discuss the constraints and compare our results to the previously 
made measurements. 


2. Data 


Eor our ana lysis we have used the data from the PCA ( jJahoda 
jet al.pOOb] ) instrument on board of the Rossi X-ray Timing Ex 
plover (RXTE) satellite. Our sample consists of three neutron 
stars: 4U 1702-429, 4U 1724-307, and SAX J1810.8-2609. 
These sources were selected because they have been known to 
exhibit PRE bursts in the hard state (i.e. at low accretion rate 
where the NS is thought to c ool passively without any external 
heating; [Kajava et al.||2014| ). They also show the most robust 
evolution of the normalization down to very low luminosities 
enabling us to use them as clear examples of bursts for which 
the cooling tail method can be applied to. These bursts are visu¬ 
ally selected, based on their evolution of normalization, to have a 
good agreement with the model used . We omit one already ana- 
lyzed long burst from 4U 1724-307 ( jSuleimanov et al.|2011a|b| ) 
as there is evidence that this particular burst might have a high 
metallicity content in the atmosphere (see Appendix |A|). 

RXTE/FCA data were analyzed with the help of heasoft 
package (version 6.16) and response matrices were generated 
using PCARSP (11.7) task of this package. All of the data were 
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Fig. 1. Reduced distributions for the black body spectral fits con¬ 
sisting of the points used in the cooling tail analysis. The dashed curve 
shows the theoretical expected distribution. Both obtained and the¬ 
oretical distribution are normalized so that the encapsulated area of the 
curves is unity. 


fitted using xspec 12.8.1 package ( Arnaud|1996 ) where the rec¬ 
ommended systematic error of 0.5% was added to the spectra 
( Jahoda et al.|2006|). We identifie d the X-ray bursts using a simi¬ 
lar method as in Galloway et al.| ( [2008] ). The time-resolved spec¬ 
tra for the bursts were extracted using an initial integration time 
of 0.25 s. In order to maintain approximately the same signal-to- 
noise ratio the integration time was doubled every time the count 
rate decreased by a factor of ^^2. The exposures were dead-time 
corrected following the approach recommended by the instru¬ 
ment teamj^The correction resulted in a roughly 10-15 per cent 
increase in the peak fiux, with the difference decreasing quickly 
as the observed count rate declined. A standard method of re¬ 
moving a 16 s spectrum taken from prior to the burst was used 


to account for the possible background emission (Kuulkers et al. 
|2002| and references therein). This standard method assumes the 
background (i.e. mainly the persistent emission) to be constant 
during the burst, even though this might not be the case. The 
changes are, however, not significant in the cooling phase (see 
Eig. 6 in jWorpel et al.||2013| ). The differences in burst charac¬ 
teristics with and without this background subtraction was also 
checked and found to be negligible at least at high burst fluxes 
(i.e. at fluxes larger than 20 per cent of the peak fiux). These 
deadtime-corrected spectra were then fitted with a blackbody 
model multiplied by an interstellar absorption model with con¬ 
stant hydrogen column density Ah (value obtained from the lit¬ 
erature, see Table [T]). The best-fit parameters are the color tem¬ 
perature Tc and the normalization constant K = (Rbb[km]/^io )^5 


^ http://heasarc.gsfc.nasa.gov/docs/xte/recipes/ 
pca_deadtime.html 
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Table 1. X-ray bursts used in the M-R analysis. 


Source 

A?h 

(10^2 cm-2) 

obsid 

Date 

(MID) 

Kxi2/K,i‘- 

4U1702-429 

YW 

50025-01-01-00 

51781.333039 

2.2 



80033-01-01-08 

52957.629763 

2.0 



80033-01-19-04 

53211.964665 

2.2 



80033-01-20-02 

53212.794286 

2.1 



80033-01-21-00 

53311.806086 

2.1 

4U 1724-307 

0.78" 

93044-06-04-00 

54526.679905 

2.3 

SAXJ1810.8-2609 

035^^ 

93044-02-04-00 

54325.894492 

3.1 


^ Ratio of the blackbody norm alizations at half-touc h down flux and at the to uc hdown (see|Kajava et al.|2014|[Foutanen et al.|2014| l, 

^ Worpel et al. ( 2013 1, ^ Kuulkers et al. ( 2003| ), ^ Natalucci et al. ( |2000[ ) 


where Z)io = /)/10 kpc. From the corresponding;^^ distributions 
(see Fig.[^ of each source we also conclude that the data is suf¬ 
ficiently well described by the blackbody model. It should also 
be noted that the theoretical atmosphere model spec tra cannot be 
perfectly fit by a (diluted) blackbody model either ( [Suleimano^ 
|et al.|20lTbl |2012| ), so in reality we do not even expect the ob- 
served distribution to be close to ideal. The bolometric fiux 
was estimated using the cflux-model in the range 0.01 - 200 
keV. All error limits were obtained using error -task in xspec. 

Some of these bursts show typical characteristics of a PRE: 
A peak in the normalization after a few seconds of the ignition. 
The evolution of the observed temperature should also show the 
characteristic double-peaked structure, arising from the cooling 
of the photosphere when it expands and the subsequent heating 
when it collapses back toward the surface due to the changing 
radiation pressure. The aforementioned signs of the expansion 
also indicate that the flux has reached (or exceeded) the Edding¬ 
ton limit during the burst. The moment when the atmosphere col¬ 
lapses back to the NS surface — i.e. the normalization reaches 
its minimum value ^td and the temperature its second peak — is 
defined as the touchdown. This also marks the beginning of the 
cooling phase where the subsequent evolution is dependent on 
the spectral state of the source, i.e. if the burst occurred in the 
hard state or in the soft state. In the hard state the normalization 
rises to a nearly constant level while the flux and the tempera¬ 
ture continue to decrease for the rest of the burst. This increase 
of normalization is due to the changing color-correction factor /c 
as we approximate the emerging spectrum with a diluted black¬ 
body model as 


stant, contrary to the theory]^ It is also crucial to notice that the 
value of the normalization in the tail of the soft-state bursts is 
different from what is observed for the hard-state b ursts. In addi¬ 
tion, t he touchdown flux can also var}0(see Fig. 1 in Kajava et al. 
|2014| ). Because of these differences, the burst selection becomes 
extremely important as our model assumptions are only valid if 
the NS surfac e alone is emitting that seems to be valid o nly in the 
hard state (see |Poutanen et al.|2014[|Kajava et al.|2Q14[ for more 
information about the soft vs. hard-state burst selection). The in¬ 
creasing emission area during the PRE phase (i.e. increase in the 
blackbody normalization K) before the touchdown is mostly re¬ 
lated to the increase of the photospheric radius. In the cooling 
tail, one believes that the evolution of K happens just because of 
varying /c with the constant actual photospheric radius equal to 
the NS radius. Therefore, the ratio of the normalizations at the 
expansion phase and the cooling tail is 


K, 




/c. 


(1 + Zx)R, 


( 3 ) 


where indices e and t refer to the expansion and tail, respectively. 
By taking /c,t ~ 1.4 in the tail and /c,e ^ 2 during the expansion 
( Pavlov et al.|1991[|Suleimanov et al.|2012| ) and demanding that 


Rq > Rt (which also means (I +Ze)^e > (1 +^t)^t for R > | ^^), 
we end up with a simple PRE condition 


^>i 

K, ^ 4 - 


( 4 ) 


Fe 


Jc 


( 1 ) 


where Be is the blackbody function and is the effective tem¬ 
perature that is connected to the observed blackbody color tem¬ 
perature as Tc = fcTQf[(l +z)“^ , where 1 -fz = (1 -IGMjRc^y^F 
is the redshift. We stress here that the decrease in the color- 
correction factor during the cooling is a feature predicted by nu- 
merous atmosphere model computations ([London et al.||1986 


[Lapidus et al.||1986} |Suleimanov et al.|[26llb[ |2012| ). Conse- 


quently, the decrease of the color correction then leads into an 
increase in the (observed) normalization because ( [Penninx et al.| 
|1989[[^n Paradijs et al.|1990| ) 


K 


- 1/4 


= /cA, 


A = (Roo[km]/Z)io) 


- 1/2 


( 2 ) 


where Roo = R(1 + z) is the apparent NS radius. On the other 
hand, for the soft-state bursts the normalization is nearly con- 


When the normalizations are equal we get Rq > 2Rt 

(note also that Zt > Ze). What is remarkable with this condition is 
that the observed “expansion”, can be less than unity (compared 
to the tail) in order for the burst to have a PRE episode. The 
PRE condition can be transformed into a requirement that the 
observed peak normalization at the expansion phase Kq must be 
larger than the normalization at the touchdown ^td as both of 
them should have similar values of the color-correction factors. 
But this is equivalent to the standard criterion that K should have 
a local minimum (at the touchdown). 


^ One possible interpretation is that in the soft state the accretion 
disk continues all the way down to the NS surface forming a spread¬ 
ing/boundary layer. A combination of emission from a partly visible 
NS and from the spreading/boundary layer itself can then create time- 
evolving spectra that appear to have almost constant color-correction 


factor (ISuleimanov & Poutanen|2006|). 

^ In the sott state the inner disk may act as a mirror re flecting part 
of the burst em ission, therefore boosting the observed flux ( [Lapidus ^ 


Sunyaev 19851. 
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Fig. 2. Bolometric flux, temperature and blackbody normalization evolution during the hard-state PRE bursts. The black line shows the estimated 
bolometric flux (left-hand y-axis) in units of 10“^erg cm"^ s"^ The blue ribbon shows the Icr limits of the blackbody normalization (outer right- 
hand y-axis) in (km/lOkpc)^. Similarly, the dashed orange line shows the color-corrected angular size (Roo/Dio)^ using the same axis. The red 
ribbon show the Icr errors for blackbody color temperature (inner right-hand y-axis) in keV. Highlighted gray area marks the region of the cooling 
tail used in the fitting procedures. 


Time-resolved spectral parameters for the bursts in our sam¬ 
ple are presented in Fig. Additionally, we show the color- 
corrected angular size with the assumption fc = 2 for the evolu¬ 
tion before the touchdown. For the /c values after the touchdown, 
we use the cooling tail model fits from Sect. to correct for the 
varying color-correction factor. Because of the new PRE criteria 
we choose to keep also the bursts that show only modest pho- 
tospheric expansion in our sample (bursts #1 and #3 from 4U 
1702-429 and burst #1 from SAX J1810.8-2609). Even though 
the expansion phase in these bursts is not very long, it is ob¬ 
vious from the Fig. that the subsequent cooling phase is still 
similar (compare, for example, the bursts #1 and #2 from 4U 
1702-429). 


3. The Bayesian cooling tail method 

For our mass and radius analysis we use the cooling t ail method 
(see [Suleimanov et al.||2011a|b| and Appendix A of |Poutanen 


|et al.||2014| ). With this method the information from the whole 
cooling track after the peak of the burst is used and the ob¬ 
served cooling is compared to the theoretical evolution predicted 
by passively cooling NS atmosphere models. To relate the ob¬ 
served data and the individual masses and radii of the NSs, we 
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use Bayesian analysis (see also Qzel & Psaltis 20151. Bayes’ 


theorem can be presented simply as (see, e.g., Grinstead & Snell 

\TWJ\ 


Pr(M|D) = 


Pr(D|M)Pr(M) 

Pr(D) 


(5) 


where Pr(Al) is the prior probability of the model At without any 
additional information from the data D, Pr(D) is the prior proba¬ 
bility of the data D, Ft(D\M) is the conditional probability of the 
data D given the model At and Fr(M\D) is the conditional prob¬ 
ability of the model M given the data D. Here the last quantity 
Pr(Al|D) is what we want to obtain and it gives us the probability 
that a given model is correct given the data. We can extend the 
Bayes theorem further by having many non-overlapping models 
At/, which exhaust the total model space At. Then the relation 
can be written as 


Pr(At/|D) = 


Pr(D|At/)Pr(At/) 

2:;Pr(D|At,-)Pr(At,-)’ 


( 6 ) 


In the cooling tail method the model space consists of four pa¬ 
rameters: mass M and radius R of the NS, hydrogen mass frac¬ 
tion X in the atmosphere, and the distance D to the source. For 
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the distance D a uniform flat prior distribution is assumed with¬ 
out any restrictions. Similarly, a uniform 2-dimensional prior 
distribution is assumed for (M, R) space. We also tak e into ac¬ 
count the causality-requirement, R > 2.S24GMIc^ { Lattimer 
|2Q12| ) and limit the mass to lie between O.SM© < M < 2.5Mo. 
For the hydrogen fraction X, a Gaussian prior distribution is used 
and is discussed in more detail later on. 

The model parameters can be combined into two new param¬ 
eters, related more closely to the color-correction curve fitting. 
The first one is the Eddington fiux 


^Edd 


GMc 


Z)2/^e(l +Z)’ 


(7) 


where Kq = 0.2(1 -r X) cm^ g ^ The second parameter is related 
to the apparent (non-color corrected) angular size 0. These pa¬ 
rameters then relate our observed fiux to the (relative) luminosity 
as F/F^dd ^ L/L^dd (where LEdd is the Eddington luminosity) 
and observed blackbody normalization to the color-correction as 
^-1/4 _ note here that it is possible to assume 

uniform priors for F^dd and A (in contrast to assuming uniform 
fiat distributio n for (M, R) space ; see Appendix [B]) as was done 
previously by |Suleimanov et al.| (|2011b| |2012|); |Poutanen et al. 

( [20T^ . 

As our actual model, we u se the recently computed hot neu¬ 
tron star atmosphere models ( [Suleimanov et al.||2012| ) that ac¬ 
count for the Klein-Nishina reduction of the electron scatter- 
ing opacity using an exact rela tivistic Compton-scattering kernel 
( Poutanen & Svensson||1996| . These models give us the color- 
correction as a function of relative luminosity, fc{£ = L/LEdd)- 
The model uncertainty is taken into account by considering a 
boxcar distribution of a width of (1 + 6) x /c (where e = 0.03) 
centered around the “real” value (see [Suleimanov et~ar] ( |2011bl 
2012|) for a discussion of model uncertainties). Compositions 


considered are a pure helium (He) and a solar composition 
of H and He with sub-solar metal abundance of Z = 0.0IZ© 
(SolAOOl). It seems that Z < O.IZ© in the surface layers of the 
NS, because in the opposite case the atmosphere model predicts 
a drop of around 20% in the /c (and correspondingly in 
at E ~ O.S/^Edd ( [Suleimanov et al.||2011b| [2012[ ), which is not 
observed. Alternatively, the absence of the drop in the low lumi¬ 
nosities might be due to an extra heating because of the accretion 
that starts again after the burst. In any case, owing to these uncer¬ 
tainties in the low luminosity regime, we neglect this area from 
the fit and consider only the regime where F > 0.2Etd, where F^d 
is the fiux at the touchdown point. 

We also relax the assumption of having a fixed hydrogen 
mass fraction Z (Z = 0.738 for SolAOOl and Z = 0 for He mod¬ 
els) and use Gaussian priors around the exact values with one 
sigma tails of 0.05. Both compositions are tested for each source 
and the physically more justified value is selected. Note also that 
this selection is easy as wrong composition gives R < 6 km 
or R > 18 km. Strictly speaking this should be taken into ac¬ 
count by using atmosphere models that are computed with the 
corresponding hydrogen fractions but the models (i.e. color- 
correction factors /c) depend so weakly on this value (as our 
one sigma limits were X = or X = 0.738 + 0.05) that it 
is possibl e to neglect the effect th at this has on the F^dd and A 
(see, e.g., [Suleimanov et al.|2012[ where the difference is rela¬ 
tively small even for Z = 0 compared to X = 1). For the M, R, 
and D this, however, has some non-negligible effect that intro¬ 
duces a small scatter of about 5 per cent to the posterior distribu¬ 
tions around the “exact” value. The uncertainty in the hydrogen 
fraction is also backed up by physical arguments because for 


hydrogen-poor companions (in the case of Z = 0, i.e. He mod¬ 
els) the evolutionary models do not exclude the possibility of the 
envelope containing some fraction of H (X < 0.1 ; [Podsiadlowski 
[et al.|200^ . The value of solar ratio of H/He is relatively accu¬ 
rately measured but here the uncertainties are possible and due to 
the possible stratification on top of the NS and/or because of the 
light ashes from the previous bursts that may stratify and accu¬ 
mulate to the surface. In the end, however, one should remember 
that the value of the hydrogen mass fraction for each model is 
still just a model assumption, made on some basis. By selecting 
a Gaussian prior around the presumed value we do weaken the 
effect that this selection has but we are unable to remove it com¬ 
pletely. If no assumption for the hydrogen mass fraction would 
be made, we could not infer the radius at all. Reassuringly, how¬ 
ever, the end results do seem to gather around similar radii which 
means that our assumed X values were close to the real values. 

In our method the data D is constructed as a set of N points 
(Fi,K7^^^) obtained from the blackbody fits, starting from the 
touchdown (/ = 1) and continuing down to 0.2Etd (i = N). The 
lower-limit here is selected so that we can maximize the avail¬ 
able data (in contrast to 1 / e F^d used in previous work) as the 
theoretical models nicely follows the data. Below the 0.2Etd limit 
the background emission can start to play too important role so 
we choose to leave it out even though some of the bursts might 
follow the model even beyond this. These data points are then 
transformed into 2-dimensional probability density distributions 
Di(F, by assuming a Gaussian measurement error model. 

Next we implicitly assume that all of the data distributions Di 
are independent of each other and also independent of the model 
assumptions and prior distributions. We then assume that the 
conditional probability of the data given the model, Pr(D|AI), 
is proportional to the product over every individual probability 
Pr(D/|M) 

Pr(D|7\4)oc Y\ Pr(D,|A<). (8) 

Each separate probability Pr(D/|AI) is assumed to be propor¬ 
tional to the path-integral evaluated through the 2-dimensional 
“data space”, (F, along the color-correction curve as 

Fr[Di\M(M,R,D,X)] 

where /c,io,hi = (1 ± ^) x /c(^) are the lower and upper limits of 
the prior boxcar distribution of the color-correction /c evaluated 
at relative luminosity £ and where e = 0.03, = J^dFEdd,^) 

is the color-correction curve in (A /c) space, the Ja¬ 

cobian transforming the path from the “model space” to the ob¬ 
served “data space” (using equations ^ and (|^), and d^- is the 
line element of The path-integral is also area-normalized (or 
length-normalized if 6 = 0) with the factor N that is defined as 
the aforementioned integral ([^ without the data D. This normal¬ 
ization takes into account the variable maximum ^ that evolves 
as a function of log g. We also note that the presented Bayesian 
path-integral formalism is related to the 2-dimensional frequen- 
tist formulation of the normalized minimum distance (see eq. (2) 
in [Poutanen et al.|201^ . 

In Bayesian interference, model parameters are determined 
using a marginal estimation where the posterior probability for a 
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Fig. 3. Left panel: Combined cooling tail in the F oc L/LEdd vs 
oc /c plane with Icr error limits presented by crosses. Similarly, 
the gray crosses show the burst evolution before the touchdown. Best-fit 
theoretieal atmosphere models are shown by the blue (SolAOOl) or red 
(He) curves. Right panel: Temperature evolution of the bursts. Black- 
body color temperature is shown for each cooling tail with black 
crosses. Red (He) or blue (SolAOOl) crosses show the color-corrected 
temperatures reff(l z)"^ The dashed lines show a powerlaw with an 
index of 4 corresponding to the F oc relation. Highlighted gray area 
marks the region of the cooling tail used in the fitting procedures. 
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model parameter pj is given by 

mpjmipj) =^f Pr[£»|M] 

X dpidp2... dpj-idpj+i... dpf^^, (10) 

where Np = 4 (corresponding to M, R, D, and X) and 

V = J Pr[£)|M] Pr[7Vl]d^M, (11) 

without the model priors that determine the integration limits. 
Then the one-dimensional function Fv[pj\D](pj) represents the 
probability that the jth parameter will take a particular value 
given the observational data D. 

The best-fit atmosphere models are presented in Fig. (left 
panel). In addition, the right panel of the figure depicts the ob¬ 
served color temperatures and the corrected effective temper¬ 
atures reff(l + z)~^ for a distant observer. Here the temperature 
is seen to follow the L oc law, i.e. showing passive cooling. 

Corresponding best-fit values of the model parameters F^dd 
and A are listed in Tablealong with the Icr and 2cr confidence 
limits of the posterior distributions. After marginalizing over the 


Fig. 4. Mass-radius constraints for the sources from the hard state 
PRE bursts. Constraints are shown by 68% (dotted line) and 95% 
(solid line) confidence level contours. The upper-left region is excluded 
by constraints from the causality and general r elativistic requirements 
( [Haensel et al.|2007l|Lattimer & Prakash|2007| ). 

radius R, mass M, and hydrogen mass fraction X we get the pos¬ 
terior distribution for the distance D too, that is also listed in Ta¬ 
ble Fig. shows the 2-dimensional mass and radius probabil¬ 
ity posterior distributions of our analysis. The obtained contours 
are arched and elongated along the curves of constant Eddington 
temperatur^ 

7’Edd.co= 1 - = 1.14x10^AF‘^/K, (12) 

WsBKel l+Z 

* Eddington temperature formulated using the and A parameters 
is not strictly constant in the M - R plane because FEdd has a log g de¬ 
pendency (because of the depend ence of fc on log g). This complication 
is introduced via the new models ( |Suleimanov et al.|2012| ) that formally 
have super-Eddington luminosities due to the Klein-Nishina reduction 
of the effective cross-section. We note that our new cooling tail formal¬ 
ism allows to take this into account as we use M and R as our parameters 
(instead of ^Edd and A). 
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Table 2. Results of the Bayesian cooling tail analysis 


Source 

Model 

T’Edd 

(10"^ erg cm"^ s"^) 

A 

([km/10 kpc]-02) 

D 

(kpc) 

4U 1702-429 

4U 1724-307 

SAXJ1810.8-2609 

He 

SolAOOl 

SolAOOl 

n (+U.U2) 

^•^^-0.02 (-0.02) 

0 (+0.04) 

'^•‘^°-0.02 (-0.03) 

^ 'yQ+0.02 (+0.03) 
(-0.03) 

^ -1 Q/^+0.001 (+0.002) 

U.iyZ -0 002 (-0.004) 

^ -1 Q A+0.002 (+0.003) 

U. 154-0 003 (-0.004) 

0 1 +0-002 (+0.004) 

U.iOV -0 002 (-0.003) 

c (+0.6) 

-^-0.4 (-0.9) 

A q+0.5 (+0.7) 
^•^-0.6 (-1.1) 

A O+0.4 (+0.6) 
^••^-0.5 (-1.0) 


Notes: Errors correspond to the 68% and 95% (in parentheses) confidence levels. 


where g is the surface gravity, (Tsb is the Stefan-Boltzmann con¬ 
stant and F _7 = /^Edd/10~^ ergcm“^ s“^ The width of the con¬ 
tours are defined by the errors in TEdd,oo that originate from the 
uncertainty in F^dd, and X. Our location on this curve is de¬ 
fined by the distance D. Because the distance is free to vary in 
our analysis we end up with the aforementioned arched posteri¬ 
ors. 


Contours from 4U 1724-307 and SAX J1810.8-2609 are 
seen to be almost identical with the radius constrained be¬ 
tween about 11-13 km (for a more strict lower mass limit 
of M > 1.1 Mo) where the largest scatter is being produced by 
the unknown distance. Both of these sources are also best ex¬ 
plained by a solar-like composition with an almost zero metal- 
licity (SolAOOl model). For 4U 1702-429 the radius is seen to 
be on the same range if mass <1.8 Mq is assumed. The model 
in this case consists of a pure helium composition. 

The choice of pure helium composition causes the Bayesian 
model to favor larger mass. This happens because F^dd ^ M/(l -r 
X), so when the hydrogen fraction increases/decreases (because 
of the possible uncertainty in the X parameter) the change can 
be balanced by increasing/decreasing M as well. With the solar 
composition of H and He our X priors are symmetrical and this 
effect is canceled out. In the case of He models the lower-limit of 
X = 0 makes the hydrogen fraction prior distribution asymmet¬ 
rical and hence the bias for larger mass is present. We note, how¬ 
ever, that in addition to this bias there is a slight preference for 
the 4U 1702-429 to favor larger logg values and hence larger 
masses. Similar effect is also present in the M ys. R posteriors 
because of our choice of flat distance prior. As F^dd ^ MjD^ we 
end up oversampling the small distance values of our flat prior 
that can then create a bias that favors smaller mass. This effect 
is visible as an increased probability density around M < 1.5 for 
4U 1724-307 and SAX J1810.8-2609. Because of these model 
biases one should be careful not to give too much emphasis on 
the specific values of the Mys.R distributions presented in Fig|^ 
but to focus more on the overall structure of the posteriors given 
by the confidence contours. 


It is also possible to set some constraints to the unknown 
distance by using the cooling tail method. From the values of 
F^dd and A, inferred from the data, we can derive the maximum 
distance where we still have M and R solutions (see Appendix A 
of |Poutanen et al.|2Q14[ for the full set of equations) as 


7) 10 — 7910,max 


1.77 X 10-2 
(1 -r X)A 2 f _7 ■ 


(13) 


On the other hand, our lower-limit of the mass prior distribution 
(Mmin = 0.8 Mq) also sets a lower-limit to the distance when 
combined with F’Edd and A. From these two constraints we are 
then able to put some limits on the distance to the NS too. 


4. EoS constraints 

The final goal of the mass and radius measurements is to con¬ 
strain the pressure-density relation of the cold dense matter. Here 
we use these new mass and radius constraints from the three NSs 
to probe the EoS by applying a Bayesian analysis to the data. 

Here the model space consists of EoS parameters Pi=i^,„^Np 
in addition to the values of neutron star masses M/^i . jvm with a 
total dimensionality of our model space SiS N = Np Nm- The 
total number of neutron stars in our sample is Nm = 3 and the 
number of EoS parameters Np depends on our initial choice of 
the model (see Sect. |4.1| ). 

The data D is now constructed as a set of Nm probability dis¬ 
tributions, Di(M,R) in the (M,R) plane obtained from the cool¬ 
ing tail posteriors presented in Section Ah of these distribu¬ 
tions are normalized to unity by computing the integral 

^^high ^^high 

dM dRni(M,R) = lVL (14) 

'7Miow 

This normalization ensures that each source is treated equally in 
the analysis. As integration limits we use the same constraints 
as in the cooling tail method analysis where Miow = 0.8 M©, 
Afhigh = 2.5 Mq, Riow = 5 km and 7?high = 18 km. 

Next we again implicitly assume that ah of the new data dis¬ 
tributions Di are independent of each other and also independent 
of the model assumptions and prior distributions. We then as¬ 
sume that the conditional probability of the data given the model, 
Pr(D|A4), is proportional to the product over the probability dis¬ 
tributions Di evaluated at some mass Mi and radius (determined 
from the model) R(Mi) so that 


Pv[D\M(pi,,„,Np, Afi,...,Ar^)] 

Yl mM,R)\M- 


Mi,R=R{Mi)- 


(15) 




For the model parameters Up -r Um we assume a uniform distri¬ 
bution (i.e. weakly informative physical priors) except for a few 
physical constraints described below. 

In order to obtain ah of the posterior probability distribu¬ 
tions for the parameters we use the publicly available bamr 
code ( |Steiner||2014a| ). The TOV solver an d data analysis rou- 
tines were obtained from the O 2 SCI library ( [Steiner|[2014b| ). To 
solve the integral the code uses the Metropolis-Hastings al¬ 
gorithm to construct a Markov chain to simulate the distribution 
Vx{D\M{pi^,„^Np^ - point, an EoS parameter pi 

and the neutron star mass Mi is generated. Corresponding radius 
curve 7?(M) (and radius Ri) is then obtained by solving the TOV 
equations. From these three masses and radii, the weight func¬ 
tion Vr{D\M] is computed and the point is either accepted or 
rejected according to the Metropolis algorithm. 
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Table 3. Prior limits for the EoS parameters 


Quantity 

Lower limit 

Upper limit 

QMC parameters 


a (MeV) 

12.5 

13.5 

a 

0.47 

0.53 

S (MeV) 

29.5 

36.1 

X (MeV) 

30.0 

70.0 

Model A parameters 


ni 

0.2 

8.0 

61 (MeV fm-3) 

150 

1600 

ni 

0.2 

8.0 

62 (MeV fm-3) 

150 

1600 


0.2 

8.0 

Model C parameters 


APi (MeV/fm^) 

0 

60 

AP 2 (MeV/fm^) 

0 

300 

AF 3 (MeV/fm^) 

0 

500 

AF 4 (MeV/fm^) 

0 

500 


4.1. EoS parameterization 


E{n) = a\—\ +b\ — 
\nol \no 


+ mn 


so that we take into account all of the possible models from Gan 
dolfi et al.]p012| ) (see Table |^. On the other hand, the parame- 


( 20061 . 


ters of the second term, b and are sensitive to the high density 
physics and control the three-nucleon interactions. Furthermore, 
in our analysis we re-parameterized b and p in terms of S and 
X. Near the nuclear saturation density hq the symmetry energy 
of the neutron matter can be obtained from ( p^ as 


S = Sino) = E{n{)) - XnucC^o) = 16 MeV + a + b. 


(17) 


where £’nuc(^o) = -16 MeV is the previously mentioned nuclear 
binding energy at the saturation density. For the pressure at the 
saturation density we obtain 


X = 3^0 


dS(n) I 

dn I 


= 3(aa + bp). 


(18) 


When building our EoS we follow the work by [Steiner et al.| 
( 2015| ) and separate our density into three effective regimes: the 
crust and regions below and above the nuclear saturation density 
fiQ. We assume nuclear binding energy of £’nuc(^o) = -16 MeV 
and saturation density of 0.16 fm“^, typical values obtained from 
[Kortelainen et al.| ( |2010| ). In mass densities these values corre¬ 
spond to about 2.7 x 10^“^ gcm“^. The nuclear symmetry en¬ 
ergy (the difference between ener^ per baryon of neutron mat¬ 
ter and that of the nuclear matter^ is denoted as <S(^b), where 
is the baryon number density and S = Sino). The pressure 
of neutron-rich matter at the saturation density hq is denoted by 
X = ^noS'ino). In addition, the nuclear masses a nd theoreti¬ 
cal models im ply a correlation between X and S ( Lattimer & 
Steiner|20141) and thus we additionally restrict these parameters 
as (9.175 - 266 MeV) < X < (14.35 - 379 MeV) enclosing 
the aforementioned constraints]^ The transition density from the 
crust to core is fixed to be at nuclear baryon density of 0.08 fm“^. 
We note that the crust model (and neither the fixed transition den¬ 
sity) has almost no effect to the resulting radii as the results are 
much more dependent on the high density behavior of our EoS. 

Near the nuclear saturation density, below the core we 
employ the “ Gandolfi-Carlson-Red dy” (GCR) quantum Monte 
Carlo model ( [Gandolfi et al.|[MT^ that takes into account the 
three-body forces between the particles in the high density mat¬ 
ter. The GCR results are accurately approximated by a two poly¬ 
trope model given in terms of the energy E for the neutron matter 
at some nucleon number density n as 


We also restrict the GCR model only up to nuclear saturation 
density as the validity of the model might not hold if a phase 
transition is present. 

Above the saturation density a set of three piecewise 
polytropes are used and re ferred to as “Model A”, similar to 
[Steiner et al.| ( [2013[ [2015[ ). In this way, when parameterizing 
the high-density EoS we introduce three continuous power laws 
defining the pressure as 


P oc e 


l + l/n 


(19) 


(16) 


with coefficients (a and b) and exponents (a and P) constrained 
by QMC calculations and where is the nucleon mass. The 
parameters of the first term, a and a, are mostly sensitive to the 
low density behavior of the EoS and are responsible of the two- 
nucleon part of the interaction. The limits of a and a are chosen 


^ Quartic terms are ignored, see [Steiner ( 

^ More accurately speaking, th e constraints originate from n uclear 
masses ([Kortelainen et al.[[2010[ ), qu antum Monte Carl o model ( |Gan-[ 
[dolfi et al.|2012[ ), c hiral interactions ([Tews et al.|2013| ), and from iso- 
baric analog states (Danielewicz & Lee 2014[ |. 


as a function of the energy density e. It has been shown that it is 
possible to model a wide range of theoretical model predictions 
with these kinds of simple polytropes with a typical rms error of 
about 4% when c ompared to the actual numerical counterparts 
( Read et al.|2009| ). In practice we can mimic theoretical models 
with up to three phase transitions because they will appear as 
successive poly tropes with different indices. Model A has 5 free 
parameters: the first transition energy 6 i and the first polytrope 
index rii, the second transition energy 62 and the second poly¬ 
trope index ^ 2 , and a third polytrope index (see Table for 
the hard limits). Additionally we, of course, require that 62 > 61 
in order to avoid double-counting of the parameter space. We 
have only 5 parameters (in contrast to 6 ) because the transition 
to the first polytrope is already fixed by the EoS at the saturation 
density. 

An alternative for the high density EoS is a piecewise linear 
model referred to as “Model C” by [Steiner et al.[ ( MT3l 2Q15| ). 
Here the low density EoS is used up to 200 MeV fm“^, after 
which four line segments are considered in the 6 vs. F plane at 
fixed energy densities of 400, 600, 1000, and 1400 MeV fm“^. 
The linear relation between the two last regimes is extrapolated 
to higher densities, if needed. The Model C has four free parame¬ 
ters: APi = P{e = 400) -F (6 = 200), AF 2 = X (6 = 600)-F(6 = 400), 
AF 3 = P{e = 1000) - P{e = 600), and AF 4 = P{e = 1400) - P{e = 
1000). These pressure changes are always relative to the value of 
the pressure at the previous fixed point, effectively showing how 
sharply the pressure changes as a function of the energy density. 
This second alternative model tends to favor strong phase tran¬ 
sitions in the core so it is interesting to study the resulting dif¬ 
ferences between it and the more smoothly evolving polytropic 
model. 

In total our EoS models have 9 (QMC -h Model A) or 8 
(QMC -r Model C) free parameters. In addition to the limits set 
to the priors (summarized in Table some combinations of the 
parameters are rejected on a physical basis. More precisely, we 
ensure that all 

1. mass-radius curves produce 2Mo NS in line with the recent 

pulsar measurements ([Antoniadis et al.[2013l[Demorest et af 

[ 20101 ), 
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Table 4. Most probable values and 68 % and 95% confidence limits for the EoS parameters. 


Quantity 

95% lower limit 

68 % lower limit Most probable value / median 

68 % upper limit 

95% upper limit 

S (MeV) 

29.6 

QMC parameters (with Model A) 

30.4 32.2 

33.3 

35.0 

X (MeV) 

32.1 

42.1 

54.9 

67.7 

69.4 

ni 

0.36 

0.45 

Model A parameters 

0.55 

0.66 

0.68 

61 (MeV fm-3) 

156 

164 

712 

865 

1020 

ni 

0.25 

0.25 

0.47 

4.80 

7.55 

62 (MeV fm-3) 

531 

794 

1190 

1510 

1560 


0.95 

0.99 

1.41 

6.80 

7.76 

S (MeV) 

29.7 

QMC parameters (with Model C) 

30.4 31.8 

33.6 

35.2 

X (MeV) 

32.0 

41.4 

54.9 

68.4 

69.4 

AFi (MeV/fm^) 

5.0 

9.9 

Model C parameters 

15 

23 

31 

AF 2 (MeV/fm^) 

59 

122 

176 

194 

195 

AP 3 (MeV/fm^) 

44 

186 

345 

386 

390 

(MeV/fm^) 

12 

26 

199 

372 

385 


Notes: Eor the X and AP 4 parameters we give the median value of the flat distribution between the Icr limits. 


2. obtained EoS are causal, i.e. dPjde > 0 and 

3. EoSs are hydrodynamically stable everywhere, i.e. they have 
an increasing pressure with increasing energy density. 

In addition, during the computations, if any of the three masses 
obtained is larger than the maximum mass for the selected EoS, 
that realization is discarded and a new one is generated. 


4.2. EoS parameter results from the Bayesian analysis 

The most probable value and the corresponding Icr and 2cr lim¬ 
its for the EoSs are summari zed in Table [4| (for the com putation 
of the confidence regions, see |Lattimer & Steiner|2Q14| ). We find 
that the posterior distributions for a and a, corresponding to the 
low-density EoS behavior (which is dominated by two-body in¬ 
teractions), are almost flat. Thus, the neutron star observations 
do not constrain these parameters, as found previously ( [Steiner 


|& Gandolfi|2012| ). We find that the derivative of the nuclear sym¬ 
metry energy X is only weakly constrained. However, we do find 
a somewhat stronger constraint on the symmetry S than what has 
been obtained previously ( [Steiner et al.|2010| ). The origin of this 
constraint is the combination of the neutron star data with the 
correlation between S and X found in quantum Monte Carlo re¬ 
sults ( [Gandolfi et al.|201^ . It is also remarkable that with both 
high-density models. Model A and Model C, the symmetry en¬ 
ergy is constrained around <S 32 MeV, that is in goo d agree- 
ment with earthly measurements (<S = 28 - 34 MeV, [Klupfel[ 
et al.[2009| ). We, however, note that the parameters obtained here 


are to be considered as “local” quantities, as they are properties 
of the EoS only at densities close to the saturation density. 

Histograms for the posterior distributions of the high-density 
parameters of the Model A are presented in Eig. The index 
of the first polytrope ni is sharply peaking around 0.5 corre¬ 
sponding to a polytropic exponent 71 = 1 l/rii ~ 3. The 
large value implies a rather stiff EoS at supranuclear densities. 
This first poly trope, corresponding to the ni index, is seen to 
continue all the way up to the first transition density at 61 
700 + 150 MeV fm“^, i.e. around 4 times the saturation energy 
density 60 . In some realizations the transition is occuring already 
at around 61 200 MeV fm“^ but these cases correspond to the 

first sharp peak seen at ^2 ~ 0.5, i.e. we practically have only 2 


polytropes spanning our energy density range. In this case, the 
role of the first polytrope is superseded by the second segment 
corresponding to index ^ 2 - In the opposite case, where all three 
polytropes span the grid the second polytrope index has values 
around 1.5 (polytropic exponent 72 ~ 1-7) with a long tail ex¬ 
tending all the way up to around 8.0 (i.e. to the upper-limit of 
our prior). What this means is that the data can not constrain 
the high-density behavior of the EoS very well. As n 2 > ni, it, 
however, indicates that some softening of the EoS is present at 
higher densities. The third polytropic index is only weakly 
constrained to be > 1 (peaking around 1.5) indicating either no 
phase transitions at all or additional softening as ^3 > ^2 > m. 


Alternatively to the rather smoothly behaving polytrope 
model, we take the piecewise linear Model C that can shows 
strong phase transitions. Instead of extending the ^2 and pa¬ 
rameter upper-limits of Model A, that would create an appar¬ 
ent bias for softer EoS, we can characterize the effect of softer 
EoSs by applying this piecewise parameterization to the data, 
too. The histograms of the posterior distributions of the obtained 
pressure differences (at the fixed transition densities e = 400, 
600, 1000, and 1400 MeVfm“^) are presented in Eig. Eor 
the first segment from 6 trans = 200 to 61 = 400 MeV the 
difference is tightly constrained around APi = 15 MeV fm“^ 
with an almost symmetric Gaussian distribution with tails of 
Icr ^ 7 MeVfm“^. This introduces a strong phase transition 
to the EoS as the pressure changes by only a little. After this 
possible transition the EoS hardens out and for AP 2 and AP^ 
(at 600 and 1000 MeV fm“^) the posterior distributions are con¬ 
siderably asymmetric toward larger pressure changes peaking at 
180 and 350 MeV fm“^, respectively. Eor the final possible line 
segment, the EoS appears almost unconstrained. The sharp cut¬ 
off at higher pressures for AP 2 , AP 3 , and AP 4 appears because 
we rule out EoSs where the speed of sound exceeds the value of 
the speed of light. The modest high-value tail for AP 4 originates 
from a small group of EoSs where the central value does not ex¬ 
ceed 1000 MeV fm“^ and hence every value of the parameter is 
allowed. 
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Fig. 5. Histograms for the posterior distributions of the high-density Model A and Model C EoS parameters as implied by the data. In the first 
panel, the inset shows a magnified view of the histogram near rii ^ 0.5. Red shading corresponds to the 68 % and blue to the 95% confidence 
regions of the parameters. Limits of the figures correspond to the hard limits set to the parameter prior distributions (see Table |^. 
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Fig. 6. Obtained EoS constraints in the M - R (left panel) and in the P - 6 plane (right panel). Upper panels correspond to the QMC + Model A 
and bottom panels to the QMC + Model C EoS. Red color indicates the probability density and black lines show the 68% (dotted) and 95% (solid) 
confidence limit contours. 


4.3. Predicted EoS 

The predicted EoS obtained from the X-ray burst data is shown 
in Fig.j^in M -R and P-e planes. Each panel in the figure dis¬ 
plays an ensemble of Id-histograms over a fixed grid in one of 
the axes (note that this is not quite the same as a 2d-histogram). 
The right panels present the ensemble of histograms of the pres¬ 
sure for each energy density. This was computed in the following 
way: for each energy density, we determined the histogram bins 
of pressure which enclose 68% and 95% of the total MC weight. 
The location of those regions for each Id-histogram are outlined 
by dotted and solid curves, respectively, and these form the con¬ 
tour lines. These Icr and 2cr contour lines give constraints on 
the pressure as a function of the energy density as implied by 
the three NS data sets (see also Table |C.l| and |C.2| ). Very high 
density behavior between the Model A and C are seen to be sim¬ 
ilar as both of the models appear rather soft in this regime. At 
these very high energy density regions it is actually the maxi¬ 
mum mass requirement that constraints the pressure evolution 
(see [Steiner et al.||2013l for more extensive discussion). Most 
striking difference occur at lower energy densities where the 


sharp phase transition in the QMC -f Model C EoS in seen to 
produce large scatter to the pressure at supranuclear densities 
(around e ^ 400 MeV fm“^). 


Similarly, the left panels of Eig. [^presents our results for the 
predicted mass-radius relations. These panels present the ensem¬ 
ble of histograms of the radius over a fixed grid in neutron star 
mass with Icr and 2 cr constraints presented with dotted and solid 


lines, respectively, similar to the right panels. See also Tables C.3 


and C.4 that summarize these contour lines as well as give the 
most probable M-R curve. The width of the contours at masses 
higher than about 1.8 M© tends to be large because the available 
NS mass and radius data in our sample generally imply smaller 
masses which in turn leads to weaker constraints. The obtained 
EoS for the QMC -f Model A has a predicted radius that is al¬ 
most constant over the whole range of viable masses. The radius 
is constrained between 11.3 - 12.8 km for M = 1.4 M© {2cr 
confidence limits). Constraints this strong are obtained because 
the combination of weak (or non-existent) phase transitions and 
the NS mass and radius measurements from the cooling tail 
method compliment each other well: Cooling tail measurements 
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Table 5. Most probable values for masses and radii for NSs constrained 
to lie on one mass versus radius curve 



M 

R 

Source 

(Mo) 

(km) 

QMC + Model A 


4U1702-429 

-j O+0.2 (+0.3) 
(-0.6) 

. . q+0.4 (+0.8) 
^^•^-0.6 (-1.1) 

4U 1724-307 

^+0.4 (+0.6) 
^•^-0.3 (-0.4) 

1 o q+0.5 (+0.8) 
^^•^-0.5 (-1.0) 

SAXJ1810.8-2609 

. .+0.4 (+0.6) 
^•^-0.4 (-0.5) 

1 o (+0.8) 

^^•^-0.5 (-1.0) 

QMC + Model C 


4U1702-429 

O+0.2 (+0.3) 
^•°-0.3 (-0.7) 

. . 4+0.5 (+1.0) 

-0.5 (-1.0) 

4U 1724-307 

a+0.6 (+0.7) 
^••^-0.3 (-0.5) 

1 1 4+0.6 (+1.0) 

-0.5 (-1.1) 

SAXJ1810.8-2609 

^+0.6 (+0.9) 
^•^-0.3 (-0.4) 

1 1 ^+0.5 (+1.0) 
^^•'^-0.6 (-1.2) 

Note: Errors correspond to the 68% and 95% (in parenthesis) 
confidence levels. 


are elongated along the constant Eddington temperature curve 
that stretches from small mass and small radius to large mass and 
large radius. On the other hand, the assumption of weak phase 
transitions in the EoS forces the radius to be almost independent 
of the mass. This assumption of constant radius then eliminates 
some of the uncertainties present in the cooling tail measure¬ 
ments (mostly due to the unknown distance) as each individual 
measurement is required to have (almost) the same radius. 

With the QMC -f Model C, on the contrary, the first phase 
transition at supranuclear densities produce slightly skewed 
mass-radius curve to compensate the cooling tail burst data that 
is elongated along the constant Eddington temperature. With 
this possible phase transition present in the EoS the mass-radius 
curve is then able to support high-mass NSs with radius of about 
R ^ 11.6 km and low-mass stars with smaller radii of around 
R ^ 11.3 km simultaneously. The phase transition also causes 
a large scatter to the radius below 1 M© as the exact location 
of the turning point where the radius starts to increase again, 
cannot be constrained from the available data. Because of these 
uncertainties originating from the possible phase transition, the 
Model C shows a much larger scatter in the predicted radii at 
small masses. The radius is constrained between 10.5 - 12.5 km 
for M = lA M© (2cr confidence limits). 


4.4. Individual nnass and radius results for the NSs 

The combination of several neutron star mass and radius mea¬ 
surements with the assumption that all neutron stars must lie on 
the same mass-radius curve puts also a significant constraints 
to the mass and radius of each object. The resulting M and R 
constraints for each object are given in Eig. [7]and summarized 
in Table In general, the QMC -f Model A EoS tends to favor 
slightly larger masses and larger radii than compared to the QMC 
-h Model C. With the polytropic Model A, the resulting mass is 
tightly constrained around M ^ 1.5 M© for 4U 1724-307 and for 
SAX J1810.8-2609. Slightly larger mass of around M 1.8 M© 
is obtained for the remaining 4U 1702-429. Resulting radii are 
constrained to be around R ^ 12.0 km for each source. With the 
piecewise linear Model C the mass is about M ^ 1.3 M© for 
the 4U 1724-307 and SAX J1810.8-2609 and, again, around 
M 1.8 M© for 4U 1702-429. The obtained radii are located 
around R ^ 11.4 km. Because of the uncertainties from the pos- 
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sible phase transition occuring in the Model C, the resulting mass 
and radius constraints for each source are also much more loose. 


5. Discussion 

In this paper we have used the cooling tail method to constrain 
the mass and radius of three NS X-ray bursters: 4U 1702-429, 
4U 1724-307, and SAX J1810.8-2609. Special care was taken 
to use only the passively cooling bursts as theoretical calcu¬ 
lations of the color-correction factor /c affecting the emerging 
spectra do not take any external heating into account. In prac¬ 
tice this means that the blackbody normalization K is r equired to 
evolve during the burst (becaus e ^ fA see, e.g., Poutanen 
et al.|2014{ [Kajava et al.|2014| ) and is indeed what we observed 


for the bursts in our sample 

Eirst we introduced a new Bayesian cooling tail method and 
assumed uniform M vs. R priors in our analysis (instead of 
uniform FEdd and A priors). By marginalizing over the M, R, 
and X prior distributions we also got distance estimates for our 
sources. One advantage here is that measurements like this are 
basically done in the X-ray band where interstellar extinction 
does not play such an important role, hence reducing possible 
model dependencies originating e.g. from selection of the in¬ 
terstellar extinction model. One should, however, note that in 
our case completely different kinds of model dependencies are 
present, related, for example, to the uncertain composition of the 
accreted material (i.e. the value of the hydrogen mass fraction) or 
to the X-ray burst selection. Unfortunately, only 4U 1724-307 
has some distance measurements that we can compare against, 
as it is located in the globular cluster Terzan 2. Distance esti¬ 
mates to this source range from 5.3 to 7.7 kpc ( [Qrtolani et al. 
|1997[ using either extinction model valid for red stars or an av¬ 
erage value from some large sample, respectively) in addition 
to the more recent measurement of 7.4 + 0.5 kpc using near- 
IR observations of red giant branch stars ( [Valenti et al.| 2010). 
We note that our distance constraints are consistent with the 
lower-limit end of the measurements as Z)max ~ 6 kpc. This 
value is, however, dependent on our selection of hydrogen mass 
fraction and should not be interpreted as a strict limit. If we 
would decrease the hydrogen mass fraction (and hence increase 
the Z)max value) our resulting radii would also rapidly increase, 
creating tension between our other measurements. Interestingly, 
if we would impose a cut around 5.3 kpc into our distance 
prior, our M vs. R results would be tightly constrained around 
R = 12.0+0.3 km and M = 1.5±0.2 M© as the new distance prior 
would remove the low-mass solutions. Eor 4U 1702-429 and 
SAX J1810.8-2609 we constrained the distance to be around 
5 .6 ^q9 kpc and 4.3^^^ kpc (2cr confidence limits), respectively. 

Mass and radius constraints of 4U 1724-307 and SAX 
J1810.8-2609 are found to be almost identical with the radius 
confined between about 11-13 km. Both of these sources are 
also best modeled by a solar-like composition with almost zero 
metallicity (SolAOOl model). The best-fit model for the third, 
remaining source 4U 1702-429 consists of pure helium. This 
implies a hydrogen-poor companion like a white dwarf, that in 
turn, implies a compact binary system in order for the accre¬ 
tion to proceed via Roche lobe overflow. Best-fit values for the 
radius of this source (with X = 0) give R 13 km at around 
M = 1.5 M©, a slightly larger value compared to the two other 
sources. 

Some physical uncertainties are also still present in the X-ray 
burst M-R measurements. Eor example, no rotation is taken into 
account in the current work. Rotation affects the emerging spec¬ 
trum owing to the various special relativistic effects (see, e.g.. 
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Fig. 7. Individual mass and radius constraints for the three neutron stars used in the analysis. Left panel shows the projected mass and the middle 
panel the projected radius histograms. Red shading corresponds to the 68% and blue to the 95% confidence regions of these parameters. Right 
panels shows the full 2-dimensional mass and radius probability distributions. Contours of 68% (dotted, black line) and 95% (solid, black line) 
confidence regions are also shown. 
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Baubock et al.|2015] ) but most importantly because the radius of 
the star increases at the equator and decreases at the poles as the 
star becomes oblate ( [Morsink et al.|2007| [Baubock ct al.|2013] ) 
It, however, also introduces two new unknown free parameters to 
the fits: the spin period and the inclination. In many cases these 
parameters are not known a priori (especially the inclination). 
In addition, the flux distribution over the surface of the NS is 
unknown. The rotation does not, however, have a considerable 
effect for the M and R constraints when the spin frequency is 
moderate (< 400 Hz). Luckily, most of the X-ray bursters de¬ 
tected seem to have a frequency around this range. In our case, 
only 4U 1702-429 has a measured spin period of 329 Hz, i.e. it 
is a slow rotator and no major corrections are expected. For the 
two other sources, no burst oscillations nor persistent millisec¬ 
ond pulsations were detected and so the spin period is unknown. 
Other uncertainties may arise from the composition that can have 
an impact on the color-correction factors ( jNattila et al.||2Q15) . 
Stratification may create an almost pure hydrogen layer or, in 
the opposite situation, strong convection might mix up the whole 
photosphere thoroughly, including the burning ashes consisting 
of heavier el ements (see e.g., [Weinberg et al.|2006t [Malone et al" 


analyzed by Suleimanov et al. ( 2011a| and 


Ozel et al. 


( |2015| l. 


Ozel et al. 


Ozel et al. 


follow the NS atmosphere model predictions, whereas the soft 
state bursts never follow them (s ee Fig. [^ and also Suleimanov 


et al.[2011a} [Poutanen et al.|2014{ [Kajava et al.|2014 for discus¬ 

sion). We therefore argue that in the soft state bursts there is an 
additional physical process (i.e. the spreading boundary layer) 
acting on the burst emission, that causes the assumptions on the 
visibility of the entire NS to break down and the value of the 
color-correction factor to be different from what is predicted by 
the passively cooling atmosphere models. One should also notice 
th e completely differ ent shape of the M - R contours obtained 


[201 1| [2014[ ). There are, however, convincing implications that 
these do not affect our current measurements as the observations 
of the normalization K do really seem to follow the theoretical 
atmosphere models with the simple solar-like (or pure He) com¬ 
positions very closely. Additional confirmation is also obtained 
from the corrected temperatures that follow the F oc 7’eff law 
implying that the values of /c used are correct. In the end, the 
distance is still our biggest source of uncertainty in the measure¬ 
ments. 

Similarly to the physical uncertainties, some technical 
sources of systematic error are present. For example data se¬ 
lection for individual bursts plays a role: fitting bursts from 
the touchdown down to only half of the touchdown flux (in¬ 
stead of the 20% used here) increases the uncertainties, as one 
would expect because we use less data, but also increases the 
radius as much as about 800 m. Also by refining the cooling 
tail method (where fc - f is used in the fitting procedures) into 
a more accurate two parameter treatment (where our model is 
^-1/4 _ (wf^)f and w is the dilution factor that was previously 
assumed to be w; f~^) the radius is seen to increase by about 
500 m (Suleimanov et al., in prep). Because both of these ef¬ 
fects act to increase our radius we can understand our current 
measurements as a lower limit to the real radius. 

The results presented here constitute the first observational 
NSM-R constraints for 4U 1702-29 and SAX J1810.8-2609 
using RXTE/FCA data. In addition, we constrained the compact¬ 
ness for the NS in 4U 1724-307 that has alre ady been previously 


m |Qzel et al. ( 2015| ) where all of the M - R points are close 
to th e region where only o ne mass-radius solutions exists (see, 
e.g., [Poutanen et al.[|2QT^ for more in depth discussion about 
this), indicating that the lower-limit for the distance is close to 
the maximum distance T>max- This also creates tension for the so¬ 
lutions to exists on higher masses as the one-solution po int lies 
on the R = 4GM/c^ line (see, e.g., Ozel & Psaltis|2015] ). In our 
case, no such a tension exists, that in turn leads in to much bigger 
M - R errors. 

Specifically for the 4U 1724- 307 we obtain R ^ 12.0 + 
0.5 km whereas Ozel et al. ( 2015t obtains ~ 12.2 km. Here the 


difference is already well inside error-limits. Suleimanov et al. 
(2011 a|b[ ) have also analyzed 4U 1724-307 using a different 
hard-state burst than what is in our sample and the resulting radii 
are considerably different (R > 14 km). There is, however, a 
possibility that the atmosphere during this burst might not con¬ 
sist of hydrogen and helium only, but is enriched with the nu¬ 
clear burning ashes. These ashes would t hen affect the meas ured 
color-correction factor considerably (see [Nattila et al.|2015[ and 
Appendix [Aj for more in-depth discussion about this burst). In 
the recent analysis of 4U 1608-52 by [Poutanen et al.[ ( [2014| ) they 
also used only the hard-state bursts from the source to constrain 
the mass and the radius. The resulting radii were, again, some¬ 
what larger at around /? > 13 km for M = 1.4 M©. We, however, 
note that 4U 1608-52 is a fast rotator with an observed spin fre¬ 
quency of 620 Hz (largest observed for an X-ray burster; [Muno 
[et al.[[2002] ). This means that for an exact and reliable M - R 
analysis, a rotation-modified cooling tail method should b e used 
(Suleimanov et al., in prep|^ see also Baubock et al.|201^ 


These three measurements were then used to create our parame¬ 
terized EoS of cold dense matter. 

In general, our EoS M - R cons traints are slightly different 


when compared to Ozel et al. (2015) results as our radii tends to 
be somewhat larger in between about 10.6 - 12.4 km for Model 
C and 11.2 - 12.7 km for Model A as compared to their mea¬ 
surement of 10.1 - 11.1 km for M = 
the parameterization of the EoS in 


1.5 Mn. Note also that 


(2015) is closer to 


our QMC -r Model A polytropic formalism. One possible cause 
for the difference might be the data selection: the PRE burst we 
analyzed in this paper occurred in the hard spectral state with a 
low persistent flux level, whereas most of the bursts analyzed in 


In the second part of the paper, we combined our separate 
cooling tail measurements to obtain a parameterized EoS of the 
cold dense matter inside NSs. By utilizing a Bayesian frame¬ 
work, we studied the constraints that are possible to set to the 
EoS parameters such as the symmetry energy S and its density 
derivative X. We also relaxed our EoS priors by studying two 
different high-density models: Model A consisting of polytropes 
and Model C with piecewise linear segments in the e - P plane. 
Here the Model C with linear segments supports strong phase 
transitions in the supranuclear densities and indeed our resulting 
EoS is seen to have large transition around e ~ 400 MeV fm“^ to 
support both large radius for large-mass and smaller radius for 
small-mass stars. Model A, on the other hand, gives much tighter 
constraints as the pressure in the core evolves rather smoothly 
and hence the resulting radius is almost independent of the mass. 
This restriction of having a constant radius then asserts tight lim¬ 
its to the obtained M vs. R results. 

Nevertheless, all of the obtained constraints seem to favor 
an EoS that produces a radius in between 10.5-12.8 km (95% 
confidence limits) for a mass of 1.4 M©. One should, however. 


(2015) occur in the soft spectral state and at higher 


persistent flux. As mentioned in the introduction, hard state X- 
ray bursts - such as the ones analyzed in this paper - tend to 
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^ Preliminary results with rotating NSs and rotation-modified cooling 
tail method show a reduction of the measured radius by about 1-2 km 
for fast rotators mostly due to the oblate emitting shape of the star and 
the varying flux distribution on the surface. 
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keep in mind that in reality the aforementioned systematic errors 
might increase these limits slightly. It is also interesting to note 
that we do not necessarily need a far better precision for the M-R 
measurements as interesting constraints can already be obtained 
from the existing observations. Encouragingly, the astrophysical 
constraints also seem to agree with nuclear physics experiments, 
theoretical studies and heavy-ion experiments of neutron matter 
( [Lattimer & Lim|2QT3] ). In addition to the EoS constraints of the 
super-dense matter we also get additional constraints for individ¬ 
ual mass and radius measurements of the single NSs. By assum¬ 
ing that all of the sources must lie on the same M vs. R curve we 
can set stronger constraints than what is possible with only the 
cooling tail method alone. Especially the results with the QMC 
-h Model A shows how the unknown distance uncertainties (i.e. 
elongated M ys. R contours along the constant Eddington tem¬ 
perature) can be reduced by combining different sources with 
this kind of joint EoS fit. In the end, the mass, for example, can 
be constrained to an accuracy of about 0.2 - 0.3 M© (Icr confi¬ 
dence level) from the original limits spanning from M = 0.8 to 
2.2 M©. 

Euture prospects include an extended study of uninformative 
(for example, Jeffrey’s) priors in M - R space, both for the EoS 
and cooling tail fit parameters. This should be done to ensure 
that we do not implicitly and unknowingly transfer information 
to the mass and radius posteriors that already have many un¬ 
certainties present due to the poor measurements. Other obvious 
prospect is the addition of all possible M-R measurements in¬ 
cluding, but not restricted to, other X-ray bursting sources, qui¬ 
escent LMXBs, thermally emitting isolated NSs, neutron star 
seismology results, pulse profiles in X-ray pulsars, moment of 
inertia and crust thickness measurements. 
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Appendix A: Two hard-state burst from 4U 
1724-307 

The cooling tail method has been previously applied to another 
hard-state burst (obsid: 10090-01-01-021) from 4U 1724-307 
( Suleimanov et al.||2011a|b| ). In this case, a rather large radius 
in excess of 13 km was obtained. Since then, there has been a 
second hard-state burst from the same source that does not, how¬ 
ever, follow the same track in the plane (this new burst 

is the one we use in our an alysis; see Fig.|A.l| ) . We also note that 
the burst data presented in |Suleimanov et aL[ ( |2011b||2012| ) was 
reduced using an older RXTE data reduction pipeline without 
deadtime correction of the PCU detectors. Since 2010, the RXTE 
team also introduced a correction to the effective area of the in¬ 
strument that changed the measured flux by up to 10%. With our 
current (and final) RXTE data reduction met hods (see Sect. [2]) 
and with the new color-correction factors from ISuleimanov et ^.1 
( |2012| ) the older burst results in a radius of around 15 km at 
1.4 M© when the cooling tail method is applied with the solar 
composition model (SolAOOl). With pure hydrogen composi¬ 
tion the radius is brought down to the ~ 14 km range. We also 
note that in [Suleimancw et al.| ( |2011b| |2012| ) the first five most 
luminous points near the touchdown are ignored in the fit. Ignor¬ 
ing these points, however, leads to smaller Eddington flux and, 
hence, even larger radius of around 16-17 km at 1.4 M© with 
our current data and new analysis methods. 

The old burst also has poor values at the low luminos¬ 
ity tail (below about < 0.5Ftd) originating from an additional 
powerlaw-like distribution of high-energy photons that cannot 
be described by the blackbody law. This, in turn, might imply an 
additional heating of the surface by infalling material from the 
hot accretion flow. The extra heating does not, however, explain 
the discrepancy between the two bursts at higher luminosities. At 
these higher luminosities, it is possible that the NS atmosphere 
does not consist of hydrogen and helium only, but is enriched 
with the nuclear burning ashes ( [Weinberg et al.|20Q6| ). The pres¬ 
ence of these ashes can then have a substantial impact on the 
color-correction factor reducing it by as much as 40% ( jNattila 


et al.|2015 


. This is also supported by the exceptionally long du¬ 
ration of the burst (~ 100 s compared to ~ 30 s for the other 
bursts in our sample) that might drive up the convection. Fur¬ 
thermore, the temperature evolution of this burst also s how s ome 
deviations from the L oc law (see right panel of Fig. jA.lj ) that 
might imply inconsistencies with pure solar composition. Be¬ 
cause to these uncertainties, we choose to leave this burst out 
from our sample and only use the newer hard-state burst from 
the source that is consistent with passive cooling. 



(10 ^ erg cm ^ s M 


Fig. A.l. Left panel: Cooling tail in the E oc L/LEdd vs oc plane 
with Icr errors presented by bl ack crosses for a hard-sta te PRE burst 
from 4U 1724-307 analyzed by jSuleimanov et al!] ( |201 Ibj ). Best-fit the¬ 
oretical atmosphere model is shown with the blue curve (SolAOOl). A 
PRE burst from the same source (used in the this paper) is also presented 
with g ray crosses and gray curve . Right panel: Temperature evolution 
of the [Suleimanov et aL] ( |2011b] ) burst. Blackbody temperature Tbb is 
shown for the cooling tail with black crosses. Blue crosses show the 
color-corrected temperatures Teff(l + z)"k The straight inclined lines 
show a powerlaw with an index of 4 corresponding to the F oc rela¬ 
tion. 
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Appendix B: Uniform F^ad and A priors 

Instead of assuming uniform M ys. R priors we can assume uni¬ 
form F^dd and A priors (and compute M and R) a s was done in 
our previous work (see, e.g., Poutanen et al.|2Q14] ). This kind of 
selection of priors, however, turns out to be problematic in many 
ways. 

First of all, we cannot truly introduce log g dependency into 
to the models as we do not know M nor R beforehand. In this 
case, we are left to iteratively selecting the best-fit surface grav¬ 
ity from the three basic values (log^ = 14.0, 14.3, and 14.6) 
instead of interpolating between the models, as is done in the 
current work. 

Secondly, there exists some issues with the transformation 
from F^dd and A parameter^ to M and R as the determinant of 
the Jacobian is 


det < / 


J F^dd^A 


i- 


M,R 


cG|l-^| 


2Z)3/2/Ces7?3/2(i 


2GMX3/2' 
Rc^ ^ 


(B.l) 


From here we can see that we end up ignoring th e so-called crit- 
ical line where R = 4GM/c^, as was shown by 


Ozel & Psaltis 


( 2Q15| ). Hence, some information of M and R is already present 


in (M, R) space owing to our choice of uniform F^dd and A pri¬ 
ors. The effect from this is visible as a splitting of our M ys. R 
posteriors into two separate regions (see Fig. 1^. For clarity, 
no scatter in X or error in /c is taken into account in this case. 
The lack of (M, R) solutions near this critical line also ends up 
affecting our final joint-fit results as can be seen from Fig. |B.2| 
Because of these aforementioned issues we choose uniform pri¬ 
ors in M and R instead, as they are our final parameters that we 
want to study. We note, however, that in the cases analyzed be¬ 
fore by [Suleimanov et al.| ( |2011b| ) and [Poutanen et al.| ( |2014| ) 
the error is not critical because the solution on M - R plane is 
constrained to be far from the R = 4GM/c^ line because of the 
distance constraints. 



8 10 12 
Radius R (km) 


Fig. B.l. Mass-radius constraints for the sources from the hard state 
PRE bursts assuming uniform FEdd and A priors. Constraints are shown 
by 68% (dotted line) and 95% (solid line) confidence level contours. 
The upper-left region is excluded b y constraints from the causality 
and general relativistic requirements ([Haensel et al.|2007| |Lattimer &| 
|Prakash|2007l . 


Our definition of the secon d A pa rameter differs from the one used 


Ozel & Psaltis 


in the work by 
Conclusions, however, remain 


|2015|) where - R^ 
the same. 
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Fig. B.2. Individual mass and radius constraints for the three neutron stars used in the analysis with uniform FEdd and A priors. Left panel shows 
the projected mass and the middle panel the projected radius histograms. Red shading corresponds to the 68% and blue to the 95% confidence 
regions of the parameters. Right panels shows the full 2-dimensional mass and radius probability distributions. Contours of 68% (dotted, black 
line) and 95% (solid, black line) confidence regions are also shown. 
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Appendix C: EoS Tables and parameter correlations 


The results for the two possible parameterized EoSs are tabu¬ 
lated here for easier access. Tables [CT] and |C.2| list the pressure 
vs. density relation for the QMC-fMo del A an d QM C-hModel 
C EoSs, respectively. Similarly, Tables [C3] and |C.4] summarize 


the resulting mass vs. radius relation for the two aforementioned 
models. Additionally, we show all the EoS parameter correla¬ 


tions between each other in Eig. C.l and C.2 


Article number, page 19 of[^ 






A&A proofs: manuscript no. nsmr 


Table C.l. Most probable values and 68% and 95% confidence limits for the pressure as a function of energy density for QMC + Model A. 


Energy density 
(MeV/fm^) 

2 cr lower limit 

Icr lower limit 

Most probable value 
(MeV/fm^) 

Icr upper limit 

2 cr upper limit 

150 

1.71 

2.15 

2.64 

3.11 

3.33 

200 

3.8 

4.74 

5.73 

6.6 

7.18 

250 

6.91 

9.1 

10.4 

12.0 

13.3 

300 

12.3 

15.3 

16.9 

19.8 

22.6 

350 

19.7 

23.6 

26.0 

30.4 

35.9 

400 

29.5 

34.0 

36.5 

44.6 

55.3 

450 

42.0 

46.7 

50.4 

63.0 

81.3 

500 

57.2 

61.9 

67.8 

86.1 

109.9 

550 

74.7 

80.7 

86.7 

113.6 

139.4 

600 

94.9 

104.0 

114.1 

143.7 

169.3 

650 

117.7 

130.6 

147.6 

173.4 

199.2 

700 

143.2 

157.7 

171.3 

202.2 

229.0 

750 

169.3 

183.6 

194.7 

230.0 

258.5 

800 

192.9 

208.6 

221.3 

258.5 

288.7 

850 

213.0 

233.3 

251.5 

288.3 

319.2 

900 

230.4 

257.0 

272.2 

318.5 

350.1 

950 

246.2 

279.1 

294.3 

348.2 

381.2 

1000 

262.3 

300.2 

323.7 

377.2 

413.9 

1050 

277.6 

320.8 

344.4 

406.3 

447.0 

1100 

293.2 

341.3 

372.6 

435.7 

481.8 

1150 

308.7 

360.8 

394.7 

465.0 

517.9 

1200 

325.1 

379.9 

414.8 

494.6 

556.1 

1250 

341.0 

399.0 

448.8 

525.0 

594.5 

1300 

357.7 

418.0 

462.0 

555.8 

634.5 

1350 

373.2 

436.7 

481.1 

586.8 

674.8 

1400 

389.2 

455.6 

505.5 

618.8 

716.4 

1450 

405.1 

474.5 

531.2 

651.7 

759.3 

1500 

421.4 

493.0 

558.2 

685.1 

804.6 


Table C.2. Most probable values and 68% and 95% confidence limits for the pressure as a function of energy density for QMC + Model C. 


Energy density 
(MeV/fm^) 

2 cr lower limit 

Icr lower limit 

Most probable value 
(MeV/fm^) 

Icr upper limit 

2 cr upper limit 

150 

1.77 

2.3 

2.97 

3.35 

3.51 

200 

3.19 

4.38 

5.46 

7.32 

9.03 

250 

3.93 

6.42 

8.94 

12.3 

15.6 

300 

4.78 

8.54 

12.6 

17.4 

22.5 

350 

5.65 

10.7 

15.1 

22.6 

29.3 

400 

11.7 

17.7 

22.7 

31.5 

38.5 

450 

35.3 

47.8 

57.7 

70.1 

78.3 

500 

53.1 

81.1 

100.8 

119.0 

127.6 

550 

70.7 

114.7 

145.1 

168.0 

177.5 

600 

92.1 

146.8 

187.3 

212.5 

224.6 

650 

122.7 

165.9 

193.1 

235.8 

259.8 

700 

156.0 

181.7 

198.8 

259.8 

300.7 

750 

183.8 

199.6 

215.0 

290.1 

341.5 

800 

200.5 

224.7 

244.4 

332.2 

386.9 

850 

212.6 

254.8 

277.7 

381.5 

436.0 

900 

223.3 

285.2 

331.8 

430.8 

485.6 

950 

234.3 

316.8 

377.0 

481.3 

535.7 

1000 

249.7 

345.3 

414.7 

523.1 

578.9 

1050 

273.4 

367.9 

435.7 

546.9 

607.1 

1100 

294.9 

389.4 

457.9 

572.5 

643.0 

1150 

312.7 

409.4 

481.1 

600.4 

683.6 

1200 

328.1 

430.1 

505.5 

632.4 

726.4 

1250 

340.8 

451.1 

531.2 

669.2 

769.6 

1300 

352.4 

471.7 

563.2 

708.8 

814.1 

1350 

363.1 

492.0 

609.2 

750.4 

859.0 

1400 

373.1 

512.4 

647.6 

793.8 

904.5 

1450 

382.6 

532.8 

680.4 

838.4 

950.5 

1500 

392.3 

555.0 

715.0 

885.7 

997.0 
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Table C.3. Most probable values and 68% and 95% confidence limits for the NS radii of fixed mass for QMC + Model A. 


Mass 

(Me) 

2 cr lower limit 

Icr lower limit 

Most probable value 
(km) 

Icr upper limit 

2 cr upper limit 

0.5 

11.60 

12.05 

12.79 

13.20 

13.73 

0.6 

11.45 

11.88 

12.19 

12.95 

13.44 

0.7 

11.37 

11.81 

12.12 

12.82 

13.25 

0.8 

11.34 

11.78 

12.12 

12.73 

13.13 

0.9 

11.33 

11.78 

12.12 

12.69 

13.03 

1.0 

11.31 

11.77 

12.04 

12.63 

12.95 

1.1 

11.32 

11.75 

12.04 

12.57 

12.90 

1.2 

11.31 

11.73 

12.04 

12.52 

12.84 

1.3 

11.30 

11.71 

12.04 

12.47 

12.80 

1.4 

11.27 

11.67 

12.04 

12.42 

12.75 

1.5 

11.23 

11.62 

11.88 

12.36 

12.71 

1.6 

11.18 

11.55 

11.80 

12.30 

12.68 

1.7 

11.10 

11.45 

11.72 

12.22 

12.63 

1.8 

10.98 

11.33 

11.64 

12.14 

12.58 

1.9 

10.81 

11.17 

11.48 

12.05 

12.52 

2.0 

10.49 

10.93 

11.40 

11.95 

12.46 

2.1 

10.65 

11.03 

11.48 

12.00 

12.52 

2.2 

10.91 

11.30 

11.56 

12.20 

12.66 


Table C.4. Most probable values and 68% and 95% confidence limits for the NS radii of fixed mass for QMC + Model C. 


Mass 

(Me) 

2 cr lower limit 

Icr lower limit 

Most probable value 
(Me) 

Icr upper limit 

2 cr upper limit 

0.5 

10.17 

11.61 

12.76 

13.63 

13.87 

0.6 

10.04 

11.00 

12.68 

13.24 

13.57 

0.7 

10.02 

10.63 

11.16 

12.96 

13.41 

0.8 

10.05 

10.38 

11.06 

12.28 

13.30 

0.9 

10.10 

10.45 

11.00 

12.07 

13.19 

1.0 

10.15 

10.54 

11.02 

11.96 

13.04 

1.1 

10.22 

10.63 

11.08 

11.90 

12.84 

1.2 

10.32 

10.72 

11.16 

11.86 

12.68 

1.3 

10.42 

10.80 

11.16 

11.85 

12.56 

1.4 

10.52 

10.88 

11.24 

11.85 

12.47 

1.5 

10.59 

10.95 

11.32 

11.85 

12.41 

1.6 

10.63 

11.01 

11.32 

11.86 

12.37 

1.7 

10.63 

11.06 

11.40 

11.88 

12.34 

1.8 

10.61 

11.09 

11.40 

11.90 

12.32 

1.9 

10.53 

11.10 

11.40 

11.93 

12.31 

2.0 

10.38 

11.06 

11.45 

11.98 

12.31 

2.1 

10.45 

11.09 

11.48 

11.99 

12.30 

2.2 

10.63 

11.13 

11.48 

11.97 

12.30 
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(MeV fm-3) (MeV fm ^) 


Fig. C.l. QMC + Model A EoS parameter correlations against each other. Red coloring shows the probability density of model parameters 
obtained from our Bayesian analysis with the cooling tail data. 
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Fig. C.2. QMC + Model C EoS parameter correlations against each other. Symbols are the same as in Eig. C.l 
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